# ------------------------------------------------------------------------------------------------
### Generate Figure 17: Estimates from 10,000 simulations randomly assigned municipalities to treatment
# ------------------------------------------------------------------------------------------------

# Load data
load("../data/results/sdid_results.Rda")
load("../data/results/simulations.Rda")

# Define filter conditions for observed ATT values
filter_conditions <- list(
  turnout = list(
    outcome = "Overall Turnout",
    party = "Overall Turnout"
  ),
  pis_support = list(
    outcome = "Party Turnout", 
    party = "Law and Justice (PiS)"
  ),
  non_pis_support = list(
    outcome = "Party Turnout",
    party = "Opposition"
  )
)

# Extract observed ATT values more efficiently
obs_att <- sapply(names(filter_conditions), function(var) {
  sdid_results %>% 
    filter(
      outcome == filter_conditions[[var]]$outcome,
      party == filter_conditions[[var]]$party,
      sum_level == "Any Level",
      sum_type == "Any Type", 
      level == "Any Level",
      type == "Any Type",
      control_sample == "notyettreated",
      no300k_sample == "Full Sample",
      noproposed_sample == "Full Sample",
      sample == "Full Sample"
    ) %>%
    pull(estimate)
})

# Compute p‐values (2-sided) ----
pvals <- simulations %>%
  group_by(outcome) %>%
  summarise(
    obs = first(obs_att[
      match(outcome, c("Overall Turnout", "Government Turnout", "Opposition Turnout"))
    ]),
    p_two = mean(abs(att) >= abs(obs)),
    .groups = "drop"
  )

# Generate plot
ggplot(simulations, aes(x = att)) +
  # Density scaled to proportion of simulations
  geom_density(
    aes(y = after_stat(count / sum(count))),
    fill = "darkgrey",
    color = "black"
  ) +
  # Empirical ATT estimate as dashed vertical line
  geom_vline(
    data = pvals,
    aes(xintercept = obs),
    linetype = "dashed",
    color = "black",
    linewidth = 0.8
  ) +
  # Empirical ATT estimate annotation
  geom_text(
    data = pvals,
    aes(x = obs, y = 0),
    label = "ATT estimated on observed data",
    angle = 270,
    vjust = -0.5,
    hjust = 1.2,
    size = 6.5
  ) +
  # Faceting and styling
  facet_wrap(~outcome, nrow = 1, scales = "free_x") +
  scale_x_continuous("sDiD coefficient value") +
  scale_y_continuous(
    "Proportion of simulations",
    labels = percent_format(accuracy = 0.1)
  ) +
  theme_bw(base_size = 23) +
  theme(
    panel.grid.major = element_line(color = "grey90"),
    panel.grid.minor = element_blank(),
    strip.text = element_text(face = "bold"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.title.y = element_text(margin = margin(r = 10))
  ) +
  labs(title = NULL)

# Save the plot
ggsave("../output/plots/figure17.pdf", width = 14, height = 7, dpi = 700)
